# Write a function that takes two sims, policy and baseline,
# and plots the impacts on production & hence CO2 emissions by location 
# (Oil vs gas, US fed/nonfed, onshore vs offshore, ROW)
calc.emissions = function(sim.p, sim.b,
                          keep.vars = c('wti','hh','brent','row.gas.delivered.price.sim',
                                        'row.oil.demand','US.oil.demand','Tot.Oil.Demand',
                                        'row.gas.demand','US.gas.demand','Tot.Gas.Demand',
                                        'row.oil.supply','US.oil.supply','Tot.Oil.Supply',
                                        'row.gas.supply','US.gas.supply','Tot.Gas.Supply',
                                        'Oil.Inventories','Gas.Inventories'),
                          oil.emis.rate = 0.43, gas.emis.rate = (117+28.55)/2204.62,
                          annualize=FALSE, share.marketed=0.884)
  # oil.emis.rates = c(Nonfederal=0.43, Federal=0.43), # tCO2/bbl
  # tCO2e/mcf, including 2.36% methane leakage at 100 year GWP=34 (Raimi email on 5/4, 10:35am))
  # gas.emis.rates = c(Nonfederal=(117+28.55)/2204.62, Federal=(117+28.55)/2204.62))
{
  
  if (any(sim.p$global.sim[, date]!=sim.b$global.sim[,date])) stop("Inconsistent dates")
  if (any(sim.p$us.sim$OilProduction_AllWells[, date]!=sim.b$us.sim$OilProduction_AllWells[,date])) stop("Inconsistent dates")
  
  well.classes.string = grep('ederal', names(sim.p$us.sim$OilProduction_AllWells), v=T)
  
  oil.emis.rate.class = rep(oil.emis.rate, length(well.classes.string))
  gas.emis.rate.class = rep(gas.emis.rate, length(well.classes.string))
  
  global.diff = data.table(date=sim.p$global.sim[,date], sim.p$global.sim[,keep.vars, with=F] - sim.b$global.sim[,keep.vars, with=F])
  # global.diff.pct = data.table(date=sim.p$global.sim[,date], sim.p$global.sim[,keep.vars, with=F]/sim.b$global.sim[,keep.vars, with=F] - 1)
  
  us.p = list(liq=sim.p$us.sim$OilProduction_AllWells, gas=sim.p$us.sim$GasProduction_AllWells)
  us.b = list(liq=sim.b$us.sim$OilProduction_AllWells, gas=sim.b$us.sim$GasProduction_AllWells)
  
  us.liq.emis.p = data.table(us.p$liq[,-'date']*365.25/12*rep.row(oil.emis.rate.class, nrow(us.p$liq)) )
  us.gas.emis.p = data.table(us.p$gas[,-'date']*365.25/12*rep.row(gas.emis.rate.class, nrow(us.p$liq)) )
  us.liq.emis.b = data.table(us.b$liq[,-'date']*365.25/12*rep.row(oil.emis.rate.class, nrow(us.b$liq)) )
  us.gas.emis.b = data.table(us.b$gas[,-'date']*365.25/12*rep.row(gas.emis.rate.class, nrow(us.b$liq)) )
  
  us.emis.p = data.table(date=us.p$liq[,date], us.liq.emis.p + us.gas.emis.p)
  us.emis.b = data.table(date=us.b$liq[,date], us.liq.emis.b + us.gas.emis.b)
  setnames(us.emis.p, 'liq_Total', 'us_emis')
  setnames(us.emis.b, 'liq_Total', 'us_emis')
  
  us.emis.diff = data.table(date=us.emis.p[,date], us.emis.p[, -'date'] - us.emis.b[, -'date'])
  # oil demand is in million barrel / day. *365.25 to million barrels per year.
  # /12 to get to million barrels per month.
  # Times tCO2/barrel to get millions of tons of CO2 per month
  # gas demand is in billion cf / day. *365.25/12 to get to bcf/month, which is a million mcf.
  # Times tCO2e/mcf to get millions of tCO2 per month.
  # Note: Tot.Gas.Demand is MARKETED production, but gas.emis.rate is per mcf of GROSS
  # production (including leakage). So to calculate emissions we scale to gross by dividing by 
  # the share marketed (88.4%), then apply the emissions rate.
  global.emiss.diff = global.diff[, .(date, global_emis=Tot.Oil.Demand*365.25/12*oil.emis.rate +
                                        Tot.Gas.Demand/share.marketed*365.25/12*gas.emis.rate)]
  
  us.liq.diff = data.table(date=us.p$liq[,date], us.p$liq[,-'date'] - us.p$liq[,-'date'])
  us.liq.diff = data.table(date=us.p$liq[,date], us.p$liq[,-'date'] - us.p$liq[,-'date'])
  us.gas.diff = data.table(date=us.p$liq[,date], us.p$gas[,-'date'] - us.p$gas[,-'date'])
  us.gas.diff = data.table(date=us.p$liq[,date], us.p$gas[,-'date'] - us.p$gas[,-'date'])
  
  
  us.emis.diff[, fed.us.emis := Oil_Federal_Offshore + Gas_Federal_Offshore +
                 Oil_Federal_Onshore + Gas_Federal_Onshore]
  
  res = merge(us.emis.diff, global.emiss.diff)
  
  res
  
  if (annualize) res = res[, lapply(.SD, sum), by=year(date)]
  
  
  
  res[, global.leakage := global_emis - fed.us.emis] # how much less did global E fall than federal?
  res[, us.leakage := us_emis - fed.us.emis] # how much less did US E fall than federal?
  res[, row.leakage := global_emis - us_emis] # leakage to ROW. So global leakage = us + row
  res[, global.leakage.pct := global.leakage/-fed.us.emis]
  
  setcolorder(res, c(ifelse(annualize,'year','date'),
                     'fed.us.emis','us_emis','global_emis',
                     'global.leakage','us.leakage','row.leakage',
                     'global.leakage.pct'))
  
  return(res)
  
}

calc.deltas = function(sb, sp, yrs = c(2030, 2040, 2050), share.marketed=0.884) {
  
  price.changes = sp$mr$root / sb$mr$root - 1
  names(price.changes) = c('Oil','Gas')
  
  sb = copy(sb$sim)
  sp = copy(sp$sim)
  
  # 2030/40/50 emisisons reductions (global)
  emis.base = sb$global.sim[, .(Global.Emis.Oil = mean(Tot.Oil.Demand*365.25*oil.emis.rate.nonfed),
                                Global.Emis.Gas = mean(Tot.Gas.Demand/share.marketed*365.25*gas.emis.rate.nonfed)),
                            by=year(date)]
  emis.base[, Global.Emis.Tot := Global.Emis.Oil + Global.Emis.Gas]
  
  # US Federal
  us.base.oil = sb$us.sim$OilProduction_AllWells[, 
                                                 .(Fed.Oil= 
                                                     mean(Oil_Federal_Onshore + 
                                                            Gas_Federal_Onshore + 
                                                            Oil_Federal_Offshore+
                                                            Gas_Federal_Offshore),
                                                   US.Oil = mean(liq_Total)),
                                                 by=year(date)]
  us.base.gas = sb$us.sim$GasProduction_AllWells[, 
                                                 .(Fed.Gas= 
                                                     mean(Oil_Federal_Onshore + 
                                                            Gas_Federal_Onshore + 
                                                            Oil_Federal_Offshore+
                                                            Gas_Federal_Offshore),
                                                   US.Gas = mean(gas_Total)),
                                                 by=year(date)]
  us.base = merge(us.base.oil, us.base.gas, by='year')
  us.base[, Fed.Emis.Oil := Fed.Oil*365.25*oil.emis.rate.fed]
  us.base[, Fed.Emis.Gas := Fed.Gas*365.25*gas.emis.rate.fed]
  us.base[, Fed.Emis.Tot := Fed.Emis.Oil + Fed.Emis.Gas]
  
  us.base[, US.Emis.Oil := Fed.Emis.Oil + (US.Oil-Fed.Oil)*365.25*oil.emis.rate.nonfed]
  us.base[, US.Emis.Gas := Fed.Emis.Gas + (US.Gas-Fed.Gas)*365.25*gas.emis.rate.nonfed]
  us.base[, US.Emis.Tot := US.Emis.Oil + US.Emis.Gas]
  
  
  base = merge(emis.base, us.base, by='year')
  
  # 2030/40/50 emisisons reductions (global)
  emis.policy = sp$global.sim[, .(Global.Emis.Oil = mean(Tot.Oil.Demand*365.25*oil.emis.rate.nonfed),
                                  Global.Emis.Gas = mean(Tot.Gas.Demand/share.marketed*365.25*gas.emis.rate.nonfed)),
                              by=year(date)]
  emis.policy[, Global.Emis.Tot := Global.Emis.Oil + Global.Emis.Gas]
  
  # US Federal
  us.policy.oil = sp$us.sim$OilProduction_AllWells[, 
                                                   .(Fed.Oil= 
                                                       mean(Oil_Federal_Onshore + 
                                                              Gas_Federal_Onshore + 
                                                              Oil_Federal_Offshore+
                                                              Gas_Federal_Offshore),
                                                     US.Oil = mean(liq_Total)),
                                                   by=year(date)]
  
  us.policy.gas = sp$us.sim$GasProduction_AllWells[, 
                                                   .(Fed.Gas= 
                                                       mean(Oil_Federal_Onshore + 
                                                              Gas_Federal_Onshore + 
                                                              Oil_Federal_Offshore+
                                                              Gas_Federal_Offshore),
                                                     US.Gas = mean(gas_Total)),
                                                   by=year(date)]
  us.policy = merge(us.policy.oil, us.policy.gas, by='year')
  us.policy[, Fed.Emis.Oil := Fed.Oil*365.25*oil.emis.rate.fed]
  us.policy[, Fed.Emis.Gas := Fed.Gas*365.25*gas.emis.rate.fed]
  us.policy[, Fed.Emis.Tot := Fed.Emis.Oil + Fed.Emis.Gas]
  
  us.policy[, US.Emis.Oil := Fed.Emis.Oil + (US.Oil-Fed.Oil)*365.25*oil.emis.rate.nonfed]
  us.policy[, US.Emis.Gas := Fed.Emis.Gas + (US.Gas-Fed.Gas)*365.25*gas.emis.rate.nonfed]
  us.policy[, US.Emis.Tot := US.Emis.Oil + US.Emis.Gas]
  
  
  policy = merge(emis.policy, us.policy, by='year')
  
  deltas = data.table(as.matrix(policy)[,-1] - as.matrix(base)[,-1])
  
  deltas = rbind(deltas, t(apply(deltas[base[,year>=2020],] , 2 , mean)))
  
  deltas[, year := c(base[, as.character(year)], 'Cum. Avg.')]
  setcolorder(deltas, 'year')
  
  return(list(p.chg = price.changes,
              deltas=deltas[year %in% c(as.character(yrs), 'Cum. Avg.'),],
              leakage = deltas[year=='Cum. Avg.', 1 - Global.Emis.Tot/Fed.Emis.Tot],
              base=base[year %in% yrs],
              policy=policy[year %in% yrs]))
  
}

# Analogous function, but calculating leakage separately by product
calc.deltas.sep.leakage = function(sb, sp, yrs = c(2030, 2040, 2050), share.marketed=0.884) {
  
  price.changes = sp$mr$root / sb$mr$root - 1
  names(price.changes) = c('Oil','Gas')
  
  sb = copy(sb$sim)
  sp = copy(sp$sim)
  
  emis.base = sb$global.sim[, .(Global.Emis.Oil = mean(Tot.Oil.Demand*365.25*oil.emis.rate.nonfed),
                                Global.Emis.Gas = mean(Tot.Gas.Demand/share.marketed*365.25*gas.emis.rate.nonfed),
                                ROW.Emis.Oil = mean(row.oil.supply*365.25*oil.emis.rate.nonfed),
                                ROW.Emis.Gas = mean(row.gas.supply/share.marketed*365.25*gas.emis.rate.nonfed)),
                            by=year(date)]
  emis.base[, Global.Emis.Tot := Global.Emis.Oil + Global.Emis.Gas]
  
  prod.base = sb$global.sim[, .(Global.Prod.Oil = mean(Tot.Oil.Supply),
                                Global.Prod.Gas = mean(Tot.Gas.Supply/share.marketed),
                                ROW.Prod.Oil = mean(row.oil.supply),
                                ROW.Prod.Gas = mean(row.gas.supply/share.marketed)),
                            by=year(date)]
  
  # US Federal
  us.base.oil = sb$us.sim$OilProduction_AllWells[, 
                                                 .(Fed.Oil= 
                                                     mean(Oil_Federal_Onshore + 
                                                            Gas_Federal_Onshore + 
                                                            Oil_Federal_Offshore+
                                                            Gas_Federal_Offshore),
                                                   US.Oil = mean(liq_Total),
                                                   Fed.Oil.On= 
                                                     mean(Oil_Federal_Onshore + 
                                                            Gas_Federal_Onshore),
                                                   Fed.Oil.Off= 
                                                     mean(Oil_Federal_Offshore + 
                                                            Gas_Federal_Offshore),
                                                   Nonfed.Oil.On= 
                                                     mean(Oil_Nonfederal_Onshore + 
                                                            Gas_Nonfederal_Onshore),
                                                   Nonfed.Oil.Off= 
                                                     mean(Oil_Nonfederal_Offshore + 
                                                            Gas_Nonfederal_Offshore)),
                                                 by=year(date)]
  us.base.gas = sb$us.sim$GasProduction_AllWells[, 
                                                 .(Fed.Gas= 
                                                     mean(Oil_Federal_Onshore + 
                                                            Gas_Federal_Onshore + 
                                                            Oil_Federal_Offshore+
                                                            Gas_Federal_Offshore),
                                                   US.Gas = mean(gas_Total),
                                                   Fed.Gas.On= 
                                                     mean(Oil_Federal_Onshore + 
                                                            Gas_Federal_Onshore),
                                                   Fed.Gas.Off= 
                                                     mean(Oil_Federal_Offshore + 
                                                            Gas_Federal_Offshore),
                                                   Nonfed.Gas.On= 
                                                     mean(Oil_Nonfederal_Onshore + 
                                                            Gas_Nonfederal_Onshore),
                                                   Nonfed.Gas.Off= 
                                                     mean(Oil_Nonfederal_Offshore + 
                                                            Gas_Nonfederal_Offshore)),
                                                 by=year(date)]
  us.base = merge(us.base.oil, us.base.gas, by='year')
  us.base[, Fed.Emis.Oil := Fed.Oil*365.25*oil.emis.rate.fed]
  us.base[, Fed.Emis.Gas := Fed.Gas*365.25*gas.emis.rate.fed]
  us.base[, Fed.Emis.Tot := Fed.Emis.Oil + Fed.Emis.Gas]
  
  us.base[, US.Emis.Oil := Fed.Emis.Oil + (US.Oil-Fed.Oil)*365.25*oil.emis.rate.nonfed]
  us.base[, US.Emis.Gas := Fed.Emis.Gas + (US.Gas-Fed.Gas)*365.25*gas.emis.rate.nonfed]
  us.base[, US.Emis.Tot := US.Emis.Oil + US.Emis.Gas]
  
  
  base = merge(emis.base, us.base, by='year')
  base = merge(base, prod.base, by='year')
  
  emis.policy = sp$global.sim[, .(Global.Emis.Oil = mean(Tot.Oil.Demand*365.25*oil.emis.rate.nonfed),
                                  Global.Emis.Gas = mean(Tot.Gas.Demand/share.marketed*365.25*gas.emis.rate.nonfed),
                                  ROW.Emis.Oil = mean(row.oil.supply*365.25*oil.emis.rate.nonfed),
                                  ROW.Emis.Gas = mean(row.gas.supply/share.marketed*365.25*gas.emis.rate.nonfed)),
                              by=year(date)]
  emis.policy[, Global.Emis.Tot := Global.Emis.Oil + Global.Emis.Gas]
  
  prod.policy = sp$global.sim[, .(Global.Prod.Oil = mean(Tot.Oil.Supply),
                                  Global.Prod.Gas = mean(Tot.Gas.Supply/share.marketed),
                                  ROW.Prod.Oil = mean(row.oil.supply),
                                  ROW.Prod.Gas = mean(row.gas.supply/share.marketed)),
                              by=year(date)]
  
  # US Federal
  us.policy.oil = sp$us.sim$OilProduction_AllWells[, 
                                                   .(Fed.Oil= 
                                                       mean(Oil_Federal_Onshore + 
                                                              Gas_Federal_Onshore + 
                                                              Oil_Federal_Offshore+
                                                              Gas_Federal_Offshore),
                                                     US.Oil = mean(liq_Total),
                                                     Fed.Oil.On= 
                                                       mean(Oil_Federal_Onshore + 
                                                              Gas_Federal_Onshore),
                                                     Fed.Oil.Off= 
                                                       mean(Oil_Federal_Offshore + 
                                                              Gas_Federal_Offshore),
                                                     Nonfed.Oil.On= 
                                                       mean(Oil_Nonfederal_Onshore + 
                                                              Gas_Nonfederal_Onshore),
                                                     Nonfed.Oil.Off= 
                                                       mean(Oil_Nonfederal_Offshore + 
                                                              Gas_Nonfederal_Offshore)),
                                                   by=year(date)]
  
  us.policy.gas = sp$us.sim$GasProduction_AllWells[, 
                                                   .(Fed.Gas= 
                                                       mean(Oil_Federal_Onshore + 
                                                              Gas_Federal_Onshore + 
                                                              Oil_Federal_Offshore+
                                                              Gas_Federal_Offshore),
                                                     US.Gas = mean(gas_Total),
                                                     Fed.Gas.On= 
                                                       mean(Oil_Federal_Onshore + 
                                                              Gas_Federal_Onshore),
                                                     Fed.Gas.Off= 
                                                       mean(Oil_Federal_Offshore + 
                                                              Gas_Federal_Offshore),
                                                     Nonfed.Gas.On= 
                                                       mean(Oil_Nonfederal_Onshore + 
                                                              Gas_Nonfederal_Onshore),
                                                     Nonfed.Gas.Off= 
                                                       mean(Oil_Nonfederal_Offshore + 
                                                              Gas_Nonfederal_Offshore)),
                                                   by=year(date)]
  us.policy = merge(us.policy.oil, us.policy.gas, by='year')
  us.policy[, Fed.Emis.Oil := Fed.Oil*365.25*oil.emis.rate.fed]
  us.policy[, Fed.Emis.Gas := Fed.Gas*365.25*gas.emis.rate.fed]
  us.policy[, Fed.Emis.Tot := Fed.Emis.Oil + Fed.Emis.Gas]
  
  us.policy[, US.Emis.Oil := Fed.Emis.Oil + (US.Oil-Fed.Oil)*365.25*oil.emis.rate.nonfed]
  us.policy[, US.Emis.Gas := Fed.Emis.Gas + (US.Gas-Fed.Gas)*365.25*gas.emis.rate.nonfed]
  us.policy[, US.Emis.Tot := US.Emis.Oil + US.Emis.Gas]
  
  policy = merge(emis.policy, us.policy, by='year')
  policy = merge(policy, prod.policy, by='year')
  
  deltas = data.table(as.matrix(policy)[,-1] - as.matrix(base)[,-1])
  
  deltas = rbind(deltas, t(apply(deltas[base[,year>=2020],] , 2 , mean)))
  
  deltas[, year := c(base[, as.character(year)], 'Cum. Avg.')]
  setcolorder(deltas, 'year')
  
  return(list(p.chg = price.changes,
              deltas=deltas[year %in% c(as.character(yrs), 'Cum. Avg.'),],
              leakage = deltas[year=='Cum. Avg.', 1 - Global.Emis.Tot/Fed.Emis.Tot],
              leakage.oil = deltas[year=='Cum. Avg.', 1 - Global.Prod.Oil/Fed.Oil],
              leakage.gas = deltas[year=='Cum. Avg.', 1 - Global.Prod.Gas/Fed.Gas],
              base=base[year %in% yrs],
              policy=policy[year %in% yrs]))
  
}

plot.sim.diffs = function(sim.p, sim.b, scen.name, date.range = NULL, cex.leg=0.85, leg.y.int=0.75) {
  
  emis = calc.emissions(sim.p$sim, sim.b$sim)
  mr.list = list(mr.p=sim.p$mr, mr.b=sim.b$mr)
  sim.p = copy(sim.p$sim)
  sim.b = copy(sim.b$sim)
  # Rename and merge together relevant variables into a single data table
  keep.vars = grep('wti$|brent$|hh$|(S|s)upply|demand', names(sim.p$global.sim), v=T)
  
  setnames(sim.p$global.sim, keep.vars, paste0(keep.vars, '.policy'))
  keep.vars.policy = c('date',paste0(keep.vars, '.policy'))
  setnames(sim.b$global.sim, keep.vars, paste0(keep.vars, '.baseline'))
  keep.vars.baseline = c('date',paste0(keep.vars, '.baseline'))
  
  # Merge together global sim results
  sim = merge(sim.p$global.sim[, keep.vars.policy, with=F],
              sim.b$global.sim[, keep.vars.baseline, with=F],
              by='date')
  setcolorder(sim, c('date', sort(names(sim)[-1])))
  setcolorder(sim, c('date','wti.baseline','wti.policy'))
  p.diff = mr.list$mr.p$root/mr.list$mr.b$root - 1 # oil and gas price changes, %
  names(p.diff) = c('wti.diff','hh.diff')
  
  
  # Merge in US-specific numbers to separate federal and nonfederal US production
  fed.cols       = grep('Federal', names(sim.p$us.sim$OilProduction_AllWells))
  nonfed.cols = grep('Nonfederal', names(sim.p$us.sim$OilProduction_AllWells))
  fed.off.cols       = grep('Federal_Offshore', names(sim.p$us.sim$OilProduction_AllWells))
  nonfed.off.cols = grep('Nonfederal_Offshore', names(sim.p$us.sim$OilProduction_AllWells))
  
  
  # Calculate federal vs nonfederal oil & g production
  # Policy case, oil production
  sim.p$us.sim$OilProduction_AllWells = cbind(sim.p$us.sim$OilProduction_AllWells,
                                              Fed.Production    = sim.p$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.cols],
                                              Nonfed.Production = sim.p$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.cols],
                                              Fed.Off.Production    = sim.p$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.off.cols],
                                              Nonfed.Off.Production = sim.p$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.off.cols])
  # Policy case, gas production
  sim.p$us.sim$GasProduction_AllWells = cbind(sim.p$us.sim$GasProduction_AllWells,
                                              Fed.Production    = sim.p$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.cols],
                                              Nonfed.Production = sim.p$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.cols],
                                              Fed.Off.Production    = sim.p$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.off.cols],
                                              Nonfed.Off.Production = sim.p$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.off.cols])
  # Base case, oil production
  sim.b$us.sim$OilProduction_AllWells = cbind(sim.b$us.sim$OilProduction_AllWells,
                                              Fed.Production    = sim.b$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.cols],
                                              Nonfed.Production = sim.b$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.cols],
                                              Fed.Off.Production    = sim.b$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.off.cols],
                                              Nonfed.Off.Production = sim.b$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.off.cols])
  # Base case, gas production
  sim.b$us.sim$GasProduction_AllWells = cbind(sim.b$us.sim$GasProduction_AllWells,
                                              Fed.Production    = sim.b$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.cols],
                                              Nonfed.Production = sim.b$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.cols],
                                              Fed.Off.Production    = sim.b$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.off.cols],
                                              Nonfed.Off.Production = sim.b$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.off.cols])
  
  
  # Merge with global results
  # liq policy
  setnames(sim.p$us.sim$OilProduction_AllWells, old=names(sim.p$us.sim$OilProduction_AllWells)[-1],
           new='liq_'%&%names(sim.p$us.sim$OilProduction_AllWells)[-1]%&%'.policy')
  # gas policy
  setnames(sim.p$us.sim$GasProduction_AllWells, old=names(sim.p$us.sim$GasProduction_AllWells)[-1],
           new='gas_'%&%names(sim.p$us.sim$GasProduction_AllWells)[-1]%&%'.policy')
  # liq baseline
  setnames(sim.b$us.sim$OilProduction_AllWells, old=names(sim.b$us.sim$OilProduction_AllWells)[-1],
           new='liq_'%&%names(sim.b$us.sim$OilProduction_AllWells)[-1]%&%'.baseline')
  # gas baseline
  setnames(sim.b$us.sim$GasProduction_AllWells, old=names(sim.b$us.sim$GasProduction_AllWells)[-1],
           new='gas_'%&%names(sim.b$us.sim$GasProduction_AllWells)[-1]%&%'.baseline')
  #
  sim = merge(sim, sim.p$us.sim$OilProduction_AllWells, by='date')
  sim = merge(sim, sim.b$us.sim$OilProduction_AllWells, by='date')
  sim = merge(sim, sim.p$us.sim$GasProduction_AllWells, by='date')
  sim = merge(sim, sim.b$us.sim$GasProduction_AllWells, by='date')
  
  if (is.null(date.range)) date.range = sim[, date]
  # browser()
  sim = sim[date %in% date.range]
  emis = emis[date %in% date.range]
  
  date.grid = seq.Date(as.Date('2000-01-01'), as.Date('2060-01-01'), by='5 years')
  pdf(output.dir%&%subfolder%&%'/Summary_graph_'%&%scen.name%&%'_'%&%today()%&%'.pdf', family='serif', width=10, height=8.5)
  par(mfrow=c(3,3))
  # First, legend indicating policy vs baseline cases.
  plot.new()
  legend('center', legend=c('Policy Case','Baseline'), col='black', lwd=2, lty=1:2, bty='n',
         y.intersp = leg.y.int, cex=2)
  
  
  # Prices
  plot(brent.policy ~ date, data=sim, type='n', col='red', ylim=c(0, 120), lwd=2,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='$/barrel (left) or $/mmbtu (right)', main='Oil & Gas Prices')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0)
  lines(brent.policy ~ date, data=sim, col='red', lwd=2, lty=1)
  lines(brent.baseline ~ date, data=sim, col='red', lwd=2, lty=2)
  lines(wti.policy ~ date, data=sim, col='black', lwd=2, lty=1)
  lines(wti.baseline ~ date, data=sim, col='black', lwd=2, lty=2)
  par(new=TRUE)
  plot(hh.policy ~ date, data=sim, type='n', col='blue', ylim=c(0, 6), lwd=2,
       xaxt='n', yaxt='n', ylab='', xlab='')
  axis(4)
  lines(hh.policy ~ date, data=sim, col='blue', lwd=2, lty=1)
  lines(hh.baseline ~ date, data=sim, col='blue', lwd=2, lty=2)
  legend('topleft', y=10, legend=c('Brent (left)','WTI (left)','Henry Hub (right)'), cex=cex.leg,
         lwd=2, lty=1, col=c('red','black','blue'), bty='n', y.intersp = leg.y.int)
  
  ## Emissions, change relative to baseline
  # Global emissions
  # Note: since emis is monthly, it is monthly CO2. Multiply by 12 to get the annual rate.
  emis.range = 1.1*range(emis[,-c('date','global.leakage.pct')]*12)
  emis = emis[date>=price.start.date]
  plot(I(global_emis*12) ~ date, data=emis, type='n', col='black', ylim=emis.range, lwd=2,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='MMTCO2e per year', main='Policy Effect on\n Global CO2e Emissions (annual rate)')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0)
  lines(I(global_emis*12) ~ date, data=emis, lwd=2, col='black', lty=1)
  lines(I(fed.us.emis*12) ~ date, data=emis, lwd=2, col=my.cols[2], lty=1)
  lines(I((us_emis - fed.us.emis)*12) ~ date, data=emis, lwd=2, col=my.cols[1], lty=1)
  lines(I((global_emis - us_emis)*12) ~ date, data=emis, lwd=2, col='purple', lty=1)
  legend('bottomleft', legend=c('Global','US Federal','US Nonfederal','ROW'),
         lty=1, lwd=2, col=c('black',my.cols[2], my.cols[1], 'purple'), bty='n',
         y.intersp = leg.y.int, cex=cex.leg)
  
  ## Oil production: Total, ROW, and US 
  # (break out US into federal onshore, federal offshore, and nonfederal)
  plot(Tot.Oil.Supply.policy ~ date, data=sim, type='n', col='black', ylim=c(0, 120), lwd=2,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='mb/d', main='Global Crude Oil Production')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0)
  lines(Tot.Oil.Supply.policy ~ date, data=sim, lwd=2, col='black', lty=1)
  lines(Tot.Oil.Supply.baseline ~ date, data=sim, lwd=2, col='black', lty=2)
  lines(row.oil.supply.policy ~ date, data=sim, lwd=2, col='purple', lty=1)
  lines(row.oil.supply.baseline ~ date, data=sim, lwd=2, col='purple', lty=2)
  lines(US.oil.supply.policy ~ date, data=sim, lwd=2, col='blue', lty=1)
  lines(US.oil.supply.baseline ~ date, data=sim, lwd=2, col='blue', lty=2)
  legend('left', legend=c('Total','ROW','US'),
         lty=1, lwd=2, col=c('black','purple','blue'), bty='n',
         y.intersp = leg.y.int, cex=cex.leg)
  
  # US production breakdown
  plot(US.oil.supply.policy ~ date, data=sim, 
       type='n', col='black', ylim=c(0, 20), lwd=2,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='mb/d', main='US Crude Oil Production')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0)
  lines(US.oil.supply.policy ~ date, data=sim, lwd=2, col='blue', lty=1)
  lines(US.oil.supply.baseline ~ date, data=sim, lwd=2, col='blue', lty=2)
  lines(liq_Nonfed.Production.policy ~ date, data=sim, lwd=2, col=my.cols[1], lty=1)
  lines(liq_Nonfed.Production.baseline ~ date, data=sim, lwd=2, col=my.cols[1], lty=2)
  lines(liq_Fed.Production.policy ~ date, data=sim, lwd=2, col=my.cols[2], lty=1)
  lines(liq_Fed.Production.baseline ~ date, data=sim, lwd=2, col=my.cols[2], lty=2)
  legend('top', legend=c('US','Nonfederal','Federal'),
         lty=1, lwd=2, col=c('blue',my.cols[1], my.cols[2]), bty='n',
         y.intersp = leg.y.int, cex=cex.leg)
  
  # Federal production breakdown
  plot(US.oil.supply.policy ~ date, data=sim, 
       type='n', col='black', ylim=c(0, 4), lwd=2,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='mb/d', main='US Federal Crude Oil Production')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0)
  lines(liq_Fed.Production.policy ~ date, data=sim, lwd=2, col=my.cols[2], lty=1)
  lines(liq_Fed.Production.baseline ~ date, data=sim, lwd=2, col=my.cols[2], lty=2)
  lines(I(liq_Fed.Production.policy-liq_Fed.Off.Production.policy) ~ date, data=sim, lwd=2, col='brown', lty=1)
  lines(I(liq_Fed.Production.baseline-liq_Fed.Off.Production.baseline) ~ date, data=sim, lwd=2, col='brown', lty=2)
  lines(liq_Fed.Off.Production.policy ~ date, data=sim, lwd=2, col=my.cols[4], lty=1)
  lines(liq_Fed.Off.Production.baseline ~ date, data=sim, lwd=2, col=my.cols[4], lty=2)
  legend('top', legend=c('Federal Total','Federal Onshore','Federal Offshore'),
         lty=1, lwd=2, col=c(my.cols[2], 'brown', my.cols[4]), bty='n',
         y.intersp = leg.y.int, cex=cex.leg)
  
  ## Gas production: Total, ROW, and US 
  # (break out US into federal onshore, federal offshore, and nonfederal)
  plot(Tot.Gas.Supply.policy ~ date, data=sim, type='n', col='black', ylim=c(0, 650), lwd=2,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='bcf/d', main='Global Gas Production')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0)
  lines(Tot.Gas.Supply.policy ~ date, data=sim, lwd=2, col='black', lty=1)
  lines(Tot.Gas.Supply.baseline ~ date, data=sim, lwd=2, col='black', lty=2)
  lines(row.gas.supply.policy ~ date, data=sim, lwd=2, col='purple', lty=1)
  lines(row.gas.supply.baseline ~ date, data=sim, lwd=2, col='purple', lty=2)
  lines(US.gas.supply.policy ~ date, data=sim, lwd=2, col='blue', lty=1)
  lines(US.gas.supply.baseline ~ date, data=sim, lwd=2, col='blue', lty=2)
  legend('top', legend=c('Total','ROW','US'),
         lty=1, lwd=2, col=c('black','purple','blue'), bty='n',
         y.intersp = leg.y.int, cex=cex.leg)
  
  # US gas production breakdown
  plot(US.gas.supply.policy ~ date, data=sim, 
       type='n', col='black', ylim=c(0, 160), lwd=2,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='bcf/d', main='US Gas Production')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0)
  lines(US.gas.supply.policy ~ date, data=sim, lwd=2, col='blue', lty=1)
  lines(US.gas.supply.baseline ~ date, data=sim, lwd=2, col='blue', lty=2)
  lines(gas_Nonfed.Production.policy ~ date, data=sim, lwd=2, col=my.cols[1], lty=1)
  lines(gas_Nonfed.Production.baseline ~ date, data=sim, lwd=2, col=my.cols[1], lty=2)
  lines(gas_Fed.Production.policy ~ date, data=sim, lwd=2, col=my.cols[2], lty=1)
  lines(gas_Fed.Production.baseline ~ date, data=sim, lwd=2, col=my.cols[2], lty=2)
  legend('top', legend=c('US','Nonfederal','Federal'),
         lty=1, lwd=2, col=c('blue',my.cols[1], my.cols[2]), bty='n',
         y.intersp = leg.y.int, cex=cex.leg)
  
  # Federal production breakdown
  plot(US.gas.supply.policy ~ date, data=sim, 
       type='n', col='black', ylim=c(0, 15), lwd=2,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='bcf/d', main='US Federal Gas Production')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0)
  lines(gas_Fed.Production.policy ~ date, data=sim, lwd=2, col=my.cols[2], lty=1)
  lines(gas_Fed.Production.baseline ~ date, data=sim, lwd=2, col=my.cols[2], lty=2)
  lines(I(gas_Fed.Production.policy-gas_Fed.Off.Production.policy) ~ date, data=sim, lwd=2, col='brown', lty=1)
  lines(I(gas_Fed.Production.baseline-gas_Fed.Off.Production.baseline) ~ date, data=sim, lwd=2, col='brown', lty=2)
  lines(gas_Fed.Off.Production.policy ~ date, data=sim, lwd=2, col=my.cols[4], lty=1)
  lines(gas_Fed.Off.Production.baseline ~ date, data=sim, lwd=2, col=my.cols[4], lty=2)
  legend('top', legend=c('Federal Total','Federal Onshore','Federal Offshore'),
         lty=1, lwd=2, col=c(my.cols[2], 'brown', my.cols[4]), bty='n',
         y.intersp = leg.y.int, cex=cex.leg)
  
  dev.off()
  
}

# Function to calculate royalty and carbon adder revenue from a sim
gov.revenue = function(sim, fed.roy.share=0.6, fed.carb.share=0.6, agg=F, share.gas.marketed=0.884,
                       carb.add.growth.rate=0.02) {
  
  # Create a dt for each FEDERAL class that has time series of 
  # royalty rates, carbon adders, and production
  sim = sim$sim$us.sim
  
  fed.prices = copy(sim$PricePaths[,'Federal',,])
  
  for (wt in c('Oil','Gas')) {
    for (o in c('Onshore','Offshore')) {
      for (u in c('Covered','Uncovered')) {
        # New oil
        prod.temp = sim$OilProduction_NewWells[,c('date',paste(wt,'Federal',o, sep='_')), with=F]
        setnames(prod.temp, c('date','new_liq'))
        
        # New oil, covered
        prod.temp = merge(prod.temp,
                          sim$OilProduction_NewWellsCovered[,c('date',paste(wt,'Federal',o, sep='_')), with=F],
                          by='date')
        setnames(prod.temp, old=paste(wt,'Federal',o, sep='_'), new='new_liq_covered')
        
        # Calculate uncovered liquids
        prod.temp[, new_liq_uncovered := new_liq - new_liq_covered]
        
        # Merge in existing oil
        prod.temp = merge(prod.temp,
                          sim$OilProduction_ExWells[,c('date',paste(wt,'Federal',o, sep='_')), with=F],
                          by='date')
        setnames(prod.temp, old=paste(wt,'Federal',o, sep='_'), new='ex_liq')
        
        # Merge in new gas
        prod.temp = merge(prod.temp,
                          sim$GasProduction_NewWells[,c('date',paste(wt,'Federal',o, sep='_')), with=F],
                          by='date')
        
        setnames(prod.temp, old=paste(wt,'Federal',o, sep='_'), new='new_gas')
        
        # Merge in new gas, covered
        prod.temp = merge(prod.temp,
                          sim$GasProduction_NewWellsCovered[,c('date',paste(wt,'Federal',o, sep='_')), with=F],
                          by='date')
        setnames(prod.temp, old=paste(wt,'Federal',o, sep='_'), new='new_gas_covered')
        
        # Calculate uncovered gas
        prod.temp[, new_gas_uncovered := new_gas - new_gas_covered]
        
        # Merge in existing gas
        prod.temp = merge(prod.temp,
                          sim$GasProduction_ExWells[,c('date',paste(wt,'Federal',o, sep='_')), with=F],
                          by='date')
        setnames(prod.temp, old=paste(wt,'Federal',o, sep='_'), new='ex_gas')
        
        fed.prices[[wt,o,u]] = merge(fed.prices[[wt,o,u]][,.(date, wti, hh, roy, carb.add.per.bbl, carb.add.per.mmbtu)],
                                     prod.temp, by='date')
        rm(prod.temp)
        
        if (u=='Uncovered') {
          fed.prices[[wt,o,u]][, relevant.oil := (ex_liq+new_liq_uncovered)]
          fed.prices[[wt,o,u]][, relevant.gas := (ex_gas+new_gas_uncovered)]
        } else {
          fed.prices[[wt,o,u]][, relevant.oil := new_liq_covered]
          fed.prices[[wt,o,u]][, relevant.gas := new_gas_covered]
        }
        
        fed.prices[[wt,o,u]][, roy.rev := 
                               (wti*relevant.oil*roy +
                                  hh*share.gas.marketed*relevant.gas*1.036*roy)*365.25/1e3] 
        
        
        fed.prices[[wt,o,u]][, carb.rev := (carb.add.per.bbl*relevant.oil +
                                              carb.add.per.mmbtu*relevant.gas*1.036)*365.25/1e3]
        
        fed.prices[[wt,o,u]][, fed.roy.rev := fed.roy.share*roy.rev]
        fed.prices[[wt,o,u]][, fed.carb.rev := fed.carb.share*carb.rev]
        fed.prices[[wt,o,u]][, tot.rev := roy.rev + carb.rev]
        fed.prices[[wt,o,u]][, fed.rev := fed.roy.rev + fed.carb.rev]
      }         
    }
    
  }
  
  roy.rev.summary = do.call(cbind, lapply(fed.prices, function(x) x[,.(roy.rev)]))
  carb.rev.summary = do.call(cbind, lapply(fed.prices, function(x) x[,.(carb.rev)]))
  tot.rev.summary = do.call(cbind, lapply(fed.prices, function(x) x[,.(tot.rev)]))
  fed.rev.summary = do.call(cbind, lapply(fed.prices, function(x) x[,.(fed.rev)]))
  my.col.names = paste(rep(dimnames(fed.prices)[[1]], 2),
                       rep(dimnames(fed.prices)[[2]], each=2),
                       rep(dimnames(fed.prices)[[3]], each=4),
                       sep="_")
  
  setnames(roy.rev.summary, my.col.names)
  setnames(carb.rev.summary, my.col.names)
  setnames(tot.rev.summary, my.col.names)
  setnames(fed.rev.summary, my.col.names)
  
  roy.rev.summary[, total.revenue:=apply(.SD, 1, sum)]
  carb.rev.summary[, total.revenue:=apply(.SD, 1, sum)]
  tot.rev.summary[, total.revenue:=apply(.SD, 1, sum)]
  fed.rev.summary[, total.revenue:=apply(.SD, 1, sum)]
  
  roy.rev.summary = cbind(fed.prices[[1]][,.(date)], roy.rev.summary)
  carb.rev.summary = cbind(fed.prices[[1]][,.(date)], carb.rev.summary)
  tot.rev.summary = cbind(fed.prices[[1]][,.(date)], tot.rev.summary)
  fed.rev.summary = cbind(fed.prices[[1]][,.(date)], fed.rev.summary)
  
  if (agg) {    
    roy.rev.summary = roy.rev.summary[, lapply(.SD, mean), by=year(date), .SDcols=names(roy.rev.summary)[-1]]
    carb.rev.summary = carb.rev.summary[, lapply(.SD, mean), by=year(date), .SDcols=names(carb.rev.summary)[-1]]
    tot.rev.summary = tot.rev.summary[, lapply(.SD, mean), by=year(date), .SDcols=names(tot.rev.summary)[-1]]
    fed.rev.summary = fed.rev.summary[, lapply(.SD, mean), by=year(date), .SDcols=names(fed.rev.summary)[-1]]
  }
  return(list(roy.rev.summary = roy.rev.summary,
              carb.rev.summary = carb.rev.summary,
              tot.rev.summary = tot.rev.summary,
              fed.rev.summary = fed.rev.summary))
  
}

gov.rev.diff = function(s1, s2, agg=F) {
  r1 = gov.revenue(s1, agg=agg, share.gas.marketed=0.884)$tot.rev.summary
  r2 = gov.revenue(s2, agg=agg, share.gas.marketed=0.884)$tot.rev.summary
  return(cbind(r1[,1], r2[,-1] - r1[,-1]))
}

extract.vars.for.table = function(d, nm) {
  out = numeric(8) # oil price, gas price, fed reduc, nonfed reduc, ROW reduc, global reduc, leakage %, revenue
  names(out) = c('Oil price', 'Gas price', 'Fed US emis', 'Nonfed US emis', 'ROW emis','Global emis','Leakage rate','Roy and CA Revenue')
  out[1:2] = d$p.chg
  out[3] = d$deltas[year=='Cum. Avg.', Fed.Emis.Tot]
  out[4] = d$deltas[year=='Cum. Avg.', US.Emis.Tot - Fed.Emis.Tot]
  out[5] = d$deltas[year=='Cum. Avg.', Global.Emis.Tot - US.Emis.Tot]
  out[6] = d$deltas[year=='Cum. Avg.', Global.Emis.Tot]
  out[7] = d$leakage
  out[8] = gov.rev.diff(sim.base, sim.list[nm])[year(date) %in% 2020:2050, mean(total.revenue)]
  
  return(out)
}


# Plot 6x2 grid of federal production over time:
plot.fed.prod = function(sim.p, sim.b, date.range = NULL, cex.leg=1.4, leg.y.int=0.75, plot.leg=F, scen.title) {
  
  emis = calc.emissions(sim.p, sim.b)
  
  sim.p = copy(sim.p)
  sim.b = copy(sim.b)
  # Rename and merge together relevant variables into a single data table
  keep.vars = grep('wti$|brent$|hh$|(S|s)upply|demand', names(sim.p$global.sim), v=T)
  
  setnames(sim.p$global.sim, keep.vars, paste0(keep.vars, '.policy'))
  keep.vars.policy = c('date',paste0(keep.vars, '.policy'))
  setnames(sim.b$global.sim, keep.vars, paste0(keep.vars, '.baseline'))
  keep.vars.baseline = c('date',paste0(keep.vars, '.baseline'))
  
  # Merge together global sim results
  sim = merge(sim.p$global.sim[, keep.vars.policy, with=F],
              sim.b$global.sim[, keep.vars.baseline, with=F],
              by='date')
  setcolorder(sim, c('date', sort(names(sim)[-1])))
  setcolorder(sim, c('date','wti.baseline','wti.policy'))
  
  # Merge in US-specific numbers to separate federal and nonfederal US production
  fed.cols       = grep('Federal', names(sim.p$us.sim$OilProduction_AllWells))
  nonfed.cols = grep('Nonfederal', names(sim.p$us.sim$OilProduction_AllWells))
  fed.off.cols       = grep('Federal_Offshore', names(sim.p$us.sim$OilProduction_AllWells))
  nonfed.off.cols = grep('Nonfederal_Offshore', names(sim.p$us.sim$OilProduction_AllWells))
  
  
  # Calculate federal vs nonfederal oil & g production
  # Policy case, oil production
  sim.p$us.sim$OilProduction_AllWells = cbind(sim.p$us.sim$OilProduction_AllWells,
                                              Fed.Production    = sim.p$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.cols],
                                              Nonfed.Production = sim.p$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.cols],
                                              Fed.Off.Production    = sim.p$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.off.cols],
                                              Nonfed.Off.Production = sim.p$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.off.cols])
  # Policy case, gas production
  sim.p$us.sim$GasProduction_AllWells = cbind(sim.p$us.sim$GasProduction_AllWells,
                                              Fed.Production    = sim.p$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.cols],
                                              Nonfed.Production = sim.p$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.cols],
                                              Fed.Off.Production    = sim.p$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.off.cols],
                                              Nonfed.Off.Production = sim.p$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.off.cols])
  # Base case, oil production
  sim.b$us.sim$OilProduction_AllWells = cbind(sim.b$us.sim$OilProduction_AllWells,
                                              Fed.Production    = sim.b$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.cols],
                                              Nonfed.Production = sim.b$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.cols],
                                              Fed.Off.Production    = sim.b$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.off.cols],
                                              Nonfed.Off.Production = sim.b$us.sim$OilProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.off.cols])
  # Base case, gas production
  sim.b$us.sim$GasProduction_AllWells = cbind(sim.b$us.sim$GasProduction_AllWells,
                                              Fed.Production    = sim.b$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.cols],
                                              Nonfed.Production = sim.b$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.cols],
                                              Fed.Off.Production    = sim.b$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=fed.off.cols],
                                              Nonfed.Off.Production = sim.b$us.sim$GasProduction_AllWells[, apply(.SD, 1, sum), .SDcols=nonfed.off.cols])
  
  
  # Merge with global results
  # liq policy
  setnames(sim.p$us.sim$OilProduction_AllWells, old=names(sim.p$us.sim$OilProduction_AllWells)[-1],
           new='liq_'%&%names(sim.p$us.sim$OilProduction_AllWells)[-1]%&%'.policy')
  # gas policy
  setnames(sim.p$us.sim$GasProduction_AllWells, old=names(sim.p$us.sim$GasProduction_AllWells)[-1],
           new='gas_'%&%names(sim.p$us.sim$GasProduction_AllWells)[-1]%&%'.policy')
  # liq baseline
  setnames(sim.b$us.sim$OilProduction_AllWells, old=names(sim.b$us.sim$OilProduction_AllWells)[-1],
           new='liq_'%&%names(sim.b$us.sim$OilProduction_AllWells)[-1]%&%'.baseline')
  # gas baseline
  setnames(sim.b$us.sim$GasProduction_AllWells, old=names(sim.b$us.sim$GasProduction_AllWells)[-1],
           new='gas_'%&%names(sim.b$us.sim$GasProduction_AllWells)[-1]%&%'.baseline')
  
  sim = merge(sim, sim.p$us.sim$OilProduction_AllWells, by='date')
  sim = merge(sim, sim.b$us.sim$OilProduction_AllWells, by='date')
  sim = merge(sim, sim.p$us.sim$GasProduction_AllWells, by='date')
  sim = merge(sim, sim.b$us.sim$GasProduction_AllWells, by='date')
  
  if (is.null(date.range)) date.range = sim[, date]
  
  sim = sim[date %in% date.range]
  emis = emis[date %in% date.range]
  
  date.grid = seq.Date(as.Date('2000-01-01'), as.Date('2060-01-01'), by='5 years')
  
  # First, scenario title
  plot.new()
  text(x=0.5,y=0.1, labels=scen.title, col='black', cex=2.2, font=2)
  if (plot.leg==T) {
    legend('topleft', legend=c('Federal Total','Federal Onshore','Federal Offshore'),
           lty=1, lwd=2, col=c(my.cols[2], 'brown', my.cols[4]), bty='n',
           y.intersp = leg.y.int, cex=cex.leg, bg='white')
    legend('topright', legend=c('Baseline','Policy Case'), cex=cex.leg,
           lwd=2, col='black',lty=2:1, bg='white', bty='n')
  }
  
  # Federal oil production breakdown
  plot(US.oil.supply.policy ~ date, data=sim, 
       type='n', col='black', ylim=c(0, 4), lwd=2,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='mb/d', main='US Federal Crude Oil Production')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0)
  lines(liq_Fed.Production.policy ~ date, data=sim, lwd=3, col=my.cols[2], lty=1)
  lines(liq_Fed.Production.baseline ~ date, data=sim, lwd=3, col=my.cols[2], lty=2)
  lines(I(liq_Fed.Production.policy-liq_Fed.Off.Production.policy) ~ date, data=sim, lwd=3, col='brown', lty=1)
  lines(I(liq_Fed.Production.baseline-liq_Fed.Off.Production.baseline) ~ date, data=sim, lwd=3, col='brown', lty=2)
  lines(liq_Fed.Off.Production.policy ~ date, data=sim, lwd=3, col=my.cols[4], lty=1)
  lines(liq_Fed.Off.Production.baseline ~ date, data=sim, lwd=3, col=my.cols[4], lty=2)
  
  # Federal gas production breakdown
  plot(US.gas.supply.policy ~ date, data=sim, 
       type='n', col='black', ylim=c(0, 15), lwd=2,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='bcf/d', main='US Federal Gas Production')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0)
  lines(gas_Fed.Production.policy ~ date, data=sim, lwd=3, col=my.cols[2], lty=1)
  lines(gas_Fed.Production.baseline ~ date, data=sim, lwd=3, col=my.cols[2], lty=2)
  lines(I(gas_Fed.Production.policy-gas_Fed.Off.Production.policy) ~ date, data=sim, lwd=3, col='brown', lty=1)
  lines(I(gas_Fed.Production.baseline-gas_Fed.Off.Production.baseline) ~ date, data=sim, lwd=3, col='brown', lty=2)
  lines(gas_Fed.Off.Production.policy ~ date, data=sim, lwd=3, col=my.cols[4], lty=1)
  lines(gas_Fed.Off.Production.baseline ~ date, data=sim, lwd=3, col=my.cols[4], lty=2)
  
  
}

plot.emis.diffs = function(sim.p, sim.b, scen.name, emis.range=c(-600,400), date.range = NULL, cex.leg=1.4, leg.y.int=0.75, plot.leg=F, scen.title) {
  
  emis = calc.emissions(sim.p, sim.b)
  
  ## Emissions, change relative to baseline
  # Global emissions
  # Note: since emis is monthly, it is monthly CO2. Multiply by 12 to get the annual rate.
  # emis.range = 1.1*range(emis[,-c('date','global.leakage.pct')]*12)
  emis = emis[date %in% date.range]
  plot(I(global_emis*12) ~ date, data=emis, type='n', col='black', ylim=emis.range, lwd=3,
       cex.lab=1.2, cex.axis=1.2, cex.main=1.2, xlab='', ylab='MMTCO2e per year', main=scen.title, cex.main=1.5)
  date.grid = seq.Date(as.Date('2000-01-01'), as.Date('2060-01-01'), by='5 years')
  abline(v=date.grid, lty=3, col='lightgray')
  abline(h=seq(-600,400, by=100), lty=3, col='lightgray')
  abline(h=0)
  lines(I(global_emis*12) ~ date, data=emis, lwd=3, col='black', lty=1)
  lines(I(fed.us.emis*12) ~ date, data=emis, lwd=3, col=my.cols[2], lty=1)
  # lines(I(us_emis*12) ~ date, data=emis, lwd=3, col='blue', lty=1)
  lines(I((us_emis - fed.us.emis)*12) ~ date, data=emis, lwd=3, col=my.cols[1], lty=1)
  lines(I((global_emis - us_emis)*12) ~ date, data=emis, lwd=3, col='purple', lty=1)
  if (plot.leg) {
    legend('bottom', legend=c('Global','US Federal','US Nonfederal','ROW'),
           lty=1, lwd=3, col=c('black',my.cols[2], my.cols[1], 'purple'), bg='white', 
           y.intersp = leg.y.int, cex=cex.leg)
  }
}


calc.deltas.for.mesh = function(sb, sp, yrs = c(2030, 2040, 2050), share.marketed=0.884) {
  
  price.changes = sp$mr$root / sb$mr$root - 1
  names(price.changes) = c('Oil','Gas')
  
  sb = copy(sb$sim)
  sp = copy(sp$sim)
  
  # 2030/40/50 emisisons reductions (global)
  emis.base = sb$global.sim[, .(Global.Emis.Oil = mean(Tot.Oil.Demand*365.25*oil.emis.rate.nonfed),
                                Global.Emis.Gas = mean(Tot.Gas.Demand/share.marketed*365.25*gas.emis.rate.nonfed)),
                            by=year(date)]
  emis.base[, Global.Emis.Tot := Global.Emis.Oil + Global.Emis.Gas]
  
  prod.base = sb$global.sim[, .(Global.Prod.Oil = mean(Tot.Oil.Supply),
                                Global.Prod.Gas = mean(Tot.Gas.Supply/share.marketed)),
                            by=year(date)]
  
  # US Federal
  us.base.oil = sb$us.sim$OilProduction_AllWells[, 
                                                 .(Fed.Oil= 
                                                     mean(Oil_Federal_Onshore + 
                                                            Gas_Federal_Onshore + 
                                                            Oil_Federal_Offshore+
                                                            Gas_Federal_Offshore),
                                                   US.Oil = mean(liq_Total)),
                                                 by=year(date)]
  us.base.gas = sb$us.sim$GasProduction_AllWells[, 
                                                 .(Fed.Gas= 
                                                     mean(Oil_Federal_Onshore + 
                                                            Gas_Federal_Onshore + 
                                                            Oil_Federal_Offshore+
                                                            Gas_Federal_Offshore),
                                                   US.Gas = mean(gas_Total)),
                                                 by=year(date)]
  us.base = merge(us.base.oil, us.base.gas, by='year')
  us.base[, Fed.Emis.Oil := Fed.Oil*365.25*oil.emis.rate.fed]
  us.base[, Fed.Emis.Gas := Fed.Gas*365.25*gas.emis.rate.fed]
  us.base[, Fed.Emis.Tot := Fed.Emis.Oil + Fed.Emis.Gas]
  
  us.base[, US.Emis.Oil := Fed.Emis.Oil + (US.Oil-Fed.Oil)*365.25*oil.emis.rate.nonfed]
  us.base[, US.Emis.Gas := Fed.Emis.Gas + (US.Gas-Fed.Gas)*365.25*gas.emis.rate.nonfed]
  us.base[, US.Emis.Tot := US.Emis.Oil + US.Emis.Gas]
  
  
  base = merge(emis.base, us.base, by='year')
  base = merge(base, prod.base, by='year')
  
  # 2030/40/50 emisisons reductions (global)
  emis.policy = sp$global.sim[, .(Global.Emis.Oil = mean(Tot.Oil.Demand*365.25*oil.emis.rate.nonfed),
                                  Global.Emis.Gas = mean(Tot.Gas.Demand/share.marketed*365.25*gas.emis.rate.nonfed)),
                              by=year(date)]
  emis.policy[, Global.Emis.Tot := Global.Emis.Oil + Global.Emis.Gas]
  
  prod.policy = sp$global.sim[, .(Global.Prod.Oil = mean(Tot.Oil.Supply),
                                  Global.Prod.Gas = mean(Tot.Gas.Supply/share.marketed)),
                              by=year(date)]
  
  # US Federal
  us.policy.oil = sp$us.sim$OilProduction_AllWells[, 
                                                   .(Fed.Oil= 
                                                       mean(Oil_Federal_Onshore + 
                                                              Gas_Federal_Onshore + 
                                                              Oil_Federal_Offshore+
                                                              Gas_Federal_Offshore),
                                                     US.Oil = mean(liq_Total)),
                                                   by=year(date)]
  
  us.policy.gas = sp$us.sim$GasProduction_AllWells[, 
                                                   .(Fed.Gas= 
                                                       mean(Oil_Federal_Onshore + 
                                                              Gas_Federal_Onshore + 
                                                              Oil_Federal_Offshore+
                                                              Gas_Federal_Offshore),
                                                     US.Gas = mean(gas_Total)),
                                                   by=year(date)]
  us.policy = merge(us.policy.oil, us.policy.gas, by='year')
  us.policy[, Fed.Emis.Oil := Fed.Oil*365.25*oil.emis.rate.fed]
  us.policy[, Fed.Emis.Gas := Fed.Gas*365.25*gas.emis.rate.fed]
  us.policy[, Fed.Emis.Tot := Fed.Emis.Oil + Fed.Emis.Gas]
  
  us.policy[, US.Emis.Oil := Fed.Emis.Oil + (US.Oil-Fed.Oil)*365.25*oil.emis.rate.nonfed]
  us.policy[, US.Emis.Gas := Fed.Emis.Gas + (US.Gas-Fed.Gas)*365.25*gas.emis.rate.nonfed]
  us.policy[, US.Emis.Tot := US.Emis.Oil + US.Emis.Gas]
  
  policy = merge(emis.policy, us.policy, by='year')
  policy = merge(policy, prod.policy, by='year')
  
  deltas = data.table(as.matrix(policy)[,-1] - as.matrix(base)[,-1])
  
  deltas = rbind(deltas, t(apply(deltas[base[,year>=2020],] , 2 , mean)))
  
  deltas[, year := c(base[, as.character(year)], 'Cum. Avg.')]
  setcolorder(deltas, 'year')
  
  return(list(p.chg = price.changes,
              deltas=deltas[year %in% c(as.character(yrs), 'Cum. Avg.'),],
              leakage = deltas[year=='Cum. Avg.', 1 - Global.Emis.Tot/Fed.Emis.Tot],
              leakage.oil = deltas[year=='Cum. Avg.', 1 - Global.Prod.Oil/Fed.Oil],
              leakage.gas = deltas[year=='Cum. Avg.', 1 - Global.Prod.Gas/Fed.Gas],
              base=base[year %in% yrs],
              policy=policy[year %in% yrs]))
  
}

extract.vars.for.table.sep.leakage = function(d, nm, base=1, sim.list.temp, sim.base.temp) {
  out = numeric(10) # oil price, gas price, fed reduc, nonfed reduc, ROW reduc, global reduc, leakage %, revenue
  names(out) = c('Oil price', 'Gas price', 'Fed US emis', 'Nonfed US emis', 'ROW emis','Global emis','Leakage rate','Roy and CA Revenue')
  out[1:2] = d$p.chg
  out[3] = d$deltas[year=='Cum. Avg.', Fed.Emis.Tot]
  out[4] = d$deltas[year=='Cum. Avg.', US.Emis.Tot - Fed.Emis.Tot]
  out[5] = d$deltas[year=='Cum. Avg.', Global.Emis.Tot - US.Emis.Tot]
  out[6] = d$deltas[year=='Cum. Avg.', Global.Emis.Tot]
  out[7] = d$leakage
  out[8] = d$leakage.oil
  out[9] = d$leakage.gas
  out[10] = gov.rev.diff(sim.base.temp, sim.list.temp, agg=F)[year(date) %in% 2020:2050, mean(total.revenue)]
  
  return(out)
}

extract.vars.for.table.for.bootstrap = function(d, nm) {
  out = numeric(8) # oil price, gas price, fed reduc, nonfed reduc, ROW reduc, global reduc, leakage %, revenue
  names(out) = c('Oil price', 'Gas price', 'Fed US emis', 'Nonfed US emis', 'ROW emis','Global emis','Leakage rate','Roy and CA Revenue')
  out[1:2] = d$p.chg
  out[3] = d$deltas[year=='Cum. Avg.', Fed.Emis.Tot]
  out[4] = d$deltas[year=='Cum. Avg.', US.Emis.Tot - Fed.Emis.Tot]
  out[5] = d$deltas[year=='Cum. Avg.', Global.Emis.Tot - US.Emis.Tot]
  out[6] = d$deltas[year=='Cum. Avg.', Global.Emis.Tot]
  out[7] = d$leakage
  out[8] = gov.rev.diff(sim.base.bs, sim.list.bs[nm])[year(date) %in% 2020:2050, mean(total.revenue)]
  
  return(out)
}
